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With the aim of describing real-time electron dynamics, we introduce an adiabatic approximation 
for the equation of motion of the one-body reduced density matrix (one-matrix) . The eigenvalues 
of the one-matrix, which represent the occupation numbers of single-particle orbitals, are obtained 
from the constrained minimization of the instantaneous ground-state energy functional rather than 
from their dynamical equations. The performance of the approximation vis-a-vis nonadiabatic 
effects is assessed in real-time simulations of a two-site Hubbard model. Due to Landau-Zener- 
type transitions, the system evolves into a nonstationary state with persistent oscillations in the 
observables. The amplitude of the oscillations displays a strongly nonmonotonic dependence on the 
strength of the electron-electron interaction and the rate of variation of the external potential. We 
interpret an associated resonance behavior in the phase of the oscillations in terms of "scattering" 
with spectator energy levels. To clarify the motivation for the minimization condition, we derive 
a sequence of energy functionals Ey , for which the corresponding sequence of minimizing one- 
matrices is asymptotic to the exact one-matrix in the adiabatic limit. 

PACS numbers: 31.15.ee,31.50.Gh,71.15.Mb 



I. INTRODUCTION 

The ability to probe and control electronic states in 
molecules and nanostructures has improved consider- 
ably in recent yearsP^ These advances have highlighted 
the need for real-time simulations of strongly-driven 
and strongly-correlated electron dynamics. Among 
the methods that have been used to simulate elec- 
tron dynamics are the time-dependent Hartree-FoclP^ 
(TDHF) and multi-configuration Hartree-FockP^P ap- 
proximations, the Kadanoff-Baym equation d 11 ! 12 ! and 
Keldysh technique^ for the nonequilibrium Green's func- 
tion, and model Hamiltonian approaches. Another fam- 
ily of methods, time-dependent density functional theory 
(TDDFT) and its extensions, describes electron dynam- 
ics in terms of single-particle densities such as the particle 
density and current density. 

Density functional theories enable a favorable compro- 
mise between accuracy and the computational accessi- 
bility of systems of interest by mapping the many-body 
Schrodinger equation, a high-dimensional linear prob- 
lem, to the single-particle Kohn-Sham equations,^ a 
lower-dimensional nonlinear problem. The Runge-Gross 
(RG) theorem^ establishes the invertibility of the map- 
ping — defined via the Schrodinger equation — from 
time-dependent local external potentials v(r, t) to time- 
dependent electron densities n(r, t), given a fixed initial 
state. The invertibility of the mapping holds also in the 
nonintcracting case, which implies that the density of 
an interacting system can be reproduced by an auxiliary 
nonintcracting system with an effective potential v s (r, t) , 
under mild restrictions on the initial stated The nonin- 
teracting system is readily described by a set of single- 
particle Schrodinger equations, the time-dependent coun- 
terparts of the Kohn-Sham equations 14 in ground-state 
DFTp^l and the density can be calculated from the re- 
sulting orbitals as n(r, t) — fa \<fii(r, t)\ 2 , where fa are 



occupation numbers. 

Electronic observables of the form A = J d 3 ra(r)n(r), 
e.g. the dipolc moment, are linear functionals of 
the density, and their time development can be ob- 
tained directly from the time-dependent density. More 
general single-particle observables of the form B = 
£ CT(T , J d 3 rd 3 r'b aa >(r,r')'ipl(r)'ip a >(r'), where *w( r > r is 
a nonlocal kernel that can contain spatial derivatives, 
are also functionals of the density if the external poten- 
tial is local, but they are generally unknown functionals 
with complicated nonlinear dependence on the density. 
In contrast, such observables are linear functionals of the 
one-body reduced density matrix (one-matrix), 

■y(x,x';t) = <*(*) | ft{x')${x)\*(t)), (1) 

where x = (r, a). The one-matrix explicitly contains 
more information than the density but less informa- 
tion than the nonequilibrium single-particle Green's func- 
tion G < (xt,x't') = i($\fo(xrt)$(xt)\*) (t < t'), as 
it is the equal-time limit of the latter, j(x, x'\t) — 
—i lim t /_ >t + G < (xt,x't'). 

In this paper, we study correlated electron dynamics 
using time-dependent reduced-density-matrix functional 
theorypSHU] which is a TDDFT-like theory in which the 
basic variable is the one-matrix instead of the density. 
The equation of motion (EOM) for the one-matrix is 

id t j(xi,x' 1 ;t) = [h (xi,t) - h (x' 1 ;t)]j(x 1 ,x' 1 ;t) 

+ 2 Jdx 2 [v c (xi,x 2 ) - v c (x' 1 ,x 2 )]T(x 1 ,X2,x' 1 ,X2;t), 

(2) 

where h (x i7 t) = pf /2 + v(rj, f), v c (x 1 ,x 2 ) = |ri -r 2 | _1 , 
J dx = I C ^ r an d 

T(x 1 ,x 2 ,x' 1 ,x' 2 ;t) 

= \ m)\ ^(4)^Vi)fe)fe) \m) (3) 
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is the two-body reduced density matrix (two-matrix). We 
use atomic units, e 2 = h = m = 1, throughout the paper. 
Equation ( 2| is the first equation of the Bogoliubov-Born- 
Green-Kirkwoo d-Yvon (BBGKY) hierarchjEiH23] f 

equa- 
tions of motion, in which the equation of motion for the 
fc-body reduced density matrix contains the (k + l)-body 
reduced density matrix. It is common to "close" the hier- 
archy at some order k by expressing the (k + l)-body re- 
duced density matrix in terms of the /c-body reduced den- 
sity matrix. In principle, the hierarchy can be closed al- 
ready at the first order, Eq. because the RG theorem 
implies 18 that the two-matrix is a universal functional of 
the one-matrix, i.e., T(t) = T([j];t), given a fixed ini- 
tial state (if a vector potential is present, the statement 
follows^ from the extension of RG in Ref. HU). This 
one-matrix functional approach has been applied within 
linear respons e theo ry to calculate frequency-dependent 
polarizabilitieiP^H and excitation energies^^H f light 
diatomic molecules. In the spirit of the adiabatic approx- 
imation to the linear response equations of TDDFTpSMl 
where the frequency-dependent exchange-correlation (xc) 
kernel f xc (uj) is replaced with its static limit f xc {^> = 0), 
these calculations employed only frequency-independent 
kernels. For two-electro n sys tems, an adiabatic approx- 
imation was constructe d 27 * 28 ! that yields the full excita- 
tion spectrum exactly, including excitations of doubly- 
excited character. This represents an advantage with 
respect to TDDFT, where excitations of doubly-excited 
character are missed if a frequency-independent xc kernel 
is used.^iHsa 

In real-time simulations based on Eq. |2]), the naive 
adiabatic approximation consists in approximating the 
two-matrix functional r([7];i) at time t by the ground- 
state functional T[y] evaluated for 7 = j(t). In other 
words, the time-dependent two-matrix is approximated 
by the adiabatic extension of the ground-state two-matrix 
to the time domain. This approximation neglects the 
memory dependence of the exact functional, i.e., the ex- 
act functional r([7];i) will generally depend on 7(2') for 
all t' < t. Memory effects in real-time dynamics have 
been studied in Refs. One of the motivations for 

taking the one-matrix as basic variable is that the univer- 
sal functionals that enter the theory might have less se- 
vere memory dependence than the functionals in TDDFT 
(for a concrete example see Ref. [2T)|) . Other extensions of 
TDDFT also hope to benefit from weaker memory depen- 
dence, notably current-density functional theoryp3ES i n 
which the basic variable is the current den sity, and time- 
dependent deformation functional theory, ^ 46 * 47 ! which 
operates with a deformation tensor. 

The naive adiabatic extension approximation has not 
yet been a pplied in the one-matrix EOM for the follow- 
ing reasonP^lIl First, consider the eigenvalue equation 
/ dx' ,r y(x,x' ']t)<f>i(x' ,t) — fi(t)4>i{x,t). The eigenfunc- 
tions and eigenvalues are called natural orbitals and occu- 
pation numbers.^ The occupation number of a natural 
orbital represents its effective occupancy in the many- 
body wave function. When the adiabatic extension ap- 



proximation is applied to any of the available ground- 
state two-matrix functionals, it yields time-independent 
occupation numbers because the available funct ionals 
have an overly-restrictive forir l 19 ! 48 -! (see Sec. II A I. This 
is disappointing because although the available function- 
als are quite accurate for ground-state properties, their 
adiabatic extension misses an important aspect of the dy- 
namics even for arbitrarily slow driving. The exact time- 
dependence of the occupation numbers in model systems 
has been inferred by solving the many-body Schrodinger 
equation numerically!^ It was found that the occupation 
numbers can indeed undergo significant changes in the 
course of the time evolution, reflecting changes in the 
degree of correlation of the many-body state. 

We introduce a simple modification of the adiabatic 
extension approximation that yields time- dependent oc- 
cupation numbers even when it is applied to the available 
ground-state two-matrix functionals. In this approxima- 
tion, which we shall refer to as the instantaneous occupa- 
tion number relaxation (IONR) approximation, the nat- 
ural orbitals satisfy a time-dependent Schrodinger equa- 
tion, while the occupation numbers are obtained "011- 
the-fly" by relaxation to the minimum of an adiabatic 
energy surface. It is expected to be accurate in the adi- 
abatic regime, i.e., when E gap T ^> 1, where E gap is the 
minimum instantaneous energy gap (between the ground 
state and first excited state) and r is the characteristic 
time scale of the external potential. To clarify the mo- 
tivation for the minimization condition and estimate the 
error in the resulting occupation numbers, we carry out 
an asymptotic analysis in the adiabatic limit r — > 00. 

We evaluate the performance of the IONR approxima- 
tion by applying it to a two-site Hubbard model. Scal- 
ing r in various external potentials of the form v(r, t/r), 
we find that it is fairly accurate even beyond the region 
of validity expected from the above adiabaticity condi- 
tion. Remarkably, it is able to describe some purely 
nonadiabatic^ effects, such as Landau-Zener-type (LZ) 
transitions.^^ \\f e also assess its robustness with re- 
spect to changes in the strength of the electron-electron 
interaction, controlled by the Hubbard parameter U. The 
effect of interactions on nonadiabatic dynamics is quite 
pronounced. Varying U across a range of values does not 
result in regular, monotonic changes in the observables 
that exhibit nonadiabatic effects. 

The paper is organized as follows. In Sec. [TTJ we in- 
troduce the IONR approximation and discuss its mo- 
tivation and validity. In Sec. |HI[ the IONR approx- 
imation and the adiabatic approximation in TDDFT 
are applied to a two-site Hubbard model. We derive 
the effective Schrodinger equation for the natural or- 
bitals (Sec. Ill A 2 1 and its linea r and semilinear ver- 
sions (Sec. |IIIA3| ). In Sec. [hTc| we carry out simula- 
tions for various time-dependent external potentials. We 
compare the IONR approximation with TDHF, adiabatic 
TDDFT, and the numerically exact solution. By vary- 
ing U and r, we study the importance of correlation for 



nonadiabatic effects. Conclusions are given in Sec. IV 
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II. ADIABATIC APPROXIMATIONS FOR THE 
ONE-MATRIX EQUATION OF MOTION 

In this section, we discuss the adiabatic extension ap- 
proximation to the one-matrix EOM and introduce a 
modification in which time-dependent occupation num- 
bers are obtained "on-the-fiy" from the constrained min- 
imization of the ground-state energy functional. To clar- 
ify the motivation for this minimization condition, wc 
derive an asymptotic sequence of energy functionals in 
which the ground-state energy functional is the zeroth 
order. The corresponding sequence of minimizing one- 
matrices is asymptotic to the exact one-matrix in the 
adiabatic limit. Using the first-order energy functional, 
we estimate the error in the IONR occupation numbers. 



A. Adiabatic extension approximation 

The adiabatic extension approximation to the one- 
matrix EOM consists in approximating r([7];t) at time 
t by the ground-state reconstruction r[7] evaluated for 
7 = j(t), the self-consistently evolving one-matrix at the 
same time. The existence of r[7], which is called a re- 
construction of the two-matrix, is implied by Gilbert's 
extension^ of the Hohenberg-Kohn theorem.^ Gilbert 
proved that the ground-state wave function is a universal 
functional of the one-matrix. Approximate reconstruc- 
tions are available from many of the existing ground- 
state energy functionals in reduced-density-matrix func- 
tional theory (RDMFT), where they are used to de- 
rive approximations to the electron-electron interaction 
energy functional W[y] = Tr(VTr) (W represents the 
Coulomb interaction operator). In RDMFT, the ana- 
log of the Hohenberg-Kohn energy functional is E v [7] — 
T[y] + V[y] + W[y], where T[y] is the kinetic energy and 
V[7] = Tr(vj) is the external potential energy. 

All of the available ground-state reconstructions have 
an overly restrictive form that couples only to the two- 
index Coulomb integrals. These include the direct Uijij 
and exchange Uijji Coulomb integrals, as well as a third 
type, Uujj, where U ijk i = (4>i4> 3 \v c \ 4>k4>i) ( in this defi- 
nition the orbitals do not contain spin factors). Recon- 
structions that couple only to such Coulomb integrals 
will be called two-index reconstructions. It will be con- 
venient to express the two-matrix in the natural orbital 
basis according to 

Fijki(t) = J dx l dx 2 dx' 1 dx' 2 4>*(x 1> t)4>j(x 2 ,t) 

x r(x 1 ,x 2 ,x / 1 ,x' 2 ;t)cf) k (x / 1 ,t)(j)i(x / 2 ,t). (4) 

If the wave function is an eigenstate of S z , all of the natu- 
ral orbitals have a pure spin state, spin-up or spin-down. 
For two-index reconstructions, the most general forms 
for the spin-parallel (o~o~) and spin-anti-parallel (era) ele- 



ments ar d 56 * 57 1 

^IjkT = F^SikSji + GffSuSjk 

rgf = l-7r 6,,*,! + Gf?8 a 8 jk + HgSijSu, (5) 

where the spin indices are made explicit. Note that the 
^ijkf 7 elements are not independent because Tf^f 7 — 
—Ffjl% a . For spin-unpolarized systems, the symmetry 
properties of the two-matrix require all i*"s and G's to 
be real symmetric and H?? to be Hermitian. In all of the 
available reconstructions, the F's, G's and Hf? are taken 
to be real-valued functions of the occupation numbers.^! 
For example, in the Miiller functional F^ 7 = fia-fja-'t 

Gff = -V7Wj^<W, ffjf = 0. Using the adiabatic 
extension of Eq. ^ with real-valued F's and G's in the 
one-matrix EOM, we obtain [see Eq. ([8| below] 

3 f / fcCT =4Im Y.Yl^r'umj 

a' ijl 

= 4Im J2 H ik Ukkii- (6) 

i 

Most of the available reconstructions have Hff = as in 
the Miiller functional. Such reconstructions will be called 
JK-only^ reconstructions because they couple only to 
the direct and exchange Coulomb integrals. For such 
reconstructions the right-hand side of Eq. Q vanishes 
identically ; 19 * 25 * 48 -^ so that the occupation numbers remain 
constant in time for any external potential v(r, t), regard- 
less of its strength and rate of change. This is clearly a 
bad approximation, as even in the adiabatic regime the 
occupation numbers can undergo large net changes. 

B. Dynamical equations for the natural orbitals 
and occupation numbers 

Expressing the one-matrix EOM, Eq. (|2j, in terms of 
the natural orbitals and occup ation numbers, leads to the 
following dynamical equations! 18 * 19 * 

id t \4>k) = (i + Weft) 1 0k) , (7) 
idtfk = 4 Im.y^^Tjjki(t) (<pk4>i \v c \ 4>i<t>j) 

ijl 

= (<Pk \u\ fa) , (8) 

where t = — V 2 /2 is the kinetic energy operator and ii e g 
and u are integral operators with kernels 

v cS (xi,x[) = v(r 1 ,t)S(r 1 - r[) 
1 

+ E f ZTi < t>3( x i> t )<t>Ux r i,t), 
jk Jk ]} 

u(xi,x\) = 2 J dx 2 [v c (xi,x 2 ) - v c (x[,X2)] 

xT(x 1 ,x 2 ,x' 1 ,x 2 ;t), (9) 
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where Ujk = (4>j \u\ 4>k) and the prime indicates that the 
sum over k is restricted to k for which //. 7^ fj. As in 
the static casej^ the components i> e gjfc for any j and k 
for which fj = fk, as well as the diagonal components 
v e s,kk, are not uniquely defined; however, this is of no 
consequence as these components have no effect on the 
dynamics of the one- matrix. Changing the diagonal com- 
ponents v e s,kk leads to a redefinition of the overall time 
dependent phase of 4>k , which cancels out in the expres- 
sion 7(x, x'\ t) =Y,kfk{t)<f ) k(x,t)<t>t(x',t). Changing the 
former components simply mixes natural orbitals 4>j and 
<fik for which fj = fk , again leaving the one-matrix invari- 
ant. Such occupationally-degenerate natural orbitals are 
defined, in the first place, only up to an arbitrary unitary 
transformation within the degenerate subspace. Symme- 
try is a common source of occupational degeneracy!^! The 
natural orbitals satisfy an effective Schrodinger equation, 
while the occupation numbers satisfy a dynamical equa- 
tion in which the kinetic energy operator and external 
potential do not appear. Equations §7$^ reveal a par- 
titioning of the operator u: only its off-diagonal part, 
which is Hermitian, contributes in Eq. ([7| , while only its 
diagonal part, which is purely imaginary, contributes in 
Eq. (||). 

Equation Q can be interpreted as the single-particle 
Schrodinger equation for the natural orbitals of an aux- 
iliary noninteracting system (a generalized KS system) 
that reproduces the time-dependent one-matrix of the 
interacting system. As in the ground-state theorypSl 
the effective potential v e g is nonlocal even though the 
given external potential is local. The occupation num- 
bers are generally fractional (0 < fa < 1). Therefore, 
the generalized KS system should be interpreted as an 
ensemble stated The occupation numbers are related 
to the ensemble weights, and their time-dependence can 
be attributed to the coupling to a fictitious environment 
through the nonhcrmitian diagonal elements of u. 



C. IONR approximation 

We propose a modification of the adiabatic exten- 
sion approximation that generates time-dependent occu- 
pation numbers even for JK-only reconstructions. The 
method consists in propagating Eq. ([7]) in the adiabatic 
extension approximation, while at each instant evaluat- 
ing v e f[ at the occupation numbers that minimize the 
ground-state energy functional E v [y] subject to the con- 
straint that the natural orbitals are equal to the time- 
dependent natural orbitals at that same instant. It can 
be interpreted as an adiabatic extension with respect to 
the orbital degrees of freedom coupled to a condition of 
instantaneous relaxation to the minimum of an adiabatic 
effective energy surface for the occupation numbers. Fun- 
damentally, the reason that the dynamical equation for 
the occupation numbers is replaced by a minimum con- 
dition, which is not a differential equation in the time 
variable, is that the adiabatic limit is a singular limit 



of the time-dependent Schrodinger equation. In this pa- 
per, we formulate the approximation for systems that 
start in the ground state at the initial time. Therefore, 
the occupation numbers relax on an effective energy sur- 
face defined by the ground-state energy functional. One 
could also consider the dynamics of a system that starts 
in another adiabatic eigenstate and determine the occu- 
pation numbers from relaxation to a minimum or station- 
ary point on the corresponding adiabatic energy surface, 
provided such a surface exists locally and obeys a lo- 
cal minimum condition or stationary condition. This is 
true, for instance, for the lowest energy state of a given 
symmetry.^ An important feature of the IONR approx- 
imation is that the condition of instantaneous relaxation 
introduces a strong temporal coherence between the or- 
bitals and occupation numbers. 

The IONR approximation is conceptually similar to 
an adiabatic approximation within linear-response time- 
dependent reduced-density-matrix functional theory de- 
veloped and tested in Rcfs. It was found that 
the occupation numbers had vanishing linear response in 
the adiabatic extension approximation. Therefore, the 
static linear response equation for the occupation num- 
bers, which gives nonzero response, was extended to finite 
uj and incorporated into the frequency-dependent linear 
response equations for the natural orbitals. This assumes 
that the occupation numbers respond instantaneously to 
the time-dependent perturbation and is equivalent to the 
result that would be obtained from the linear response of 
the IONR approximation. 

D. Minimum principles and the motivation for the 
IONR approximation 

The IONR approximation can be motivated by an 
asymptotic analysis of the many-body Schrodinger equa- 
tion in the adiabatic limit r — > 00. On the basis of such 
an analysis, we identify a sequence of "ground-state" en- 
ergy functionals E^ n '[y], each of which satisfies a local 
minimum principle at each instant of time for sufficiently 

large r. The instantaneous minimization of E^ [7] gives 
an approximation y( n > with error of order <r~'™ +1 ) to 
the exact time-dependent one-matrix. In the IONR ap- 
proximation, the occupation numbers are calculated from 
the constrained minimization of the Hohenberg-Kohn- 
like functional E v [7] , which is the zeroth-order member of 
the sequence. An explicit comparison of the zeroth-order 
and first-order energy functionals affords us a means of 
estimating the error in the IONR occupation numbers. 

Consider a system with Hamiltonian of the form H = 
H(t/r) that starts in the ground state at t = —00. Fur- 
ther, suppose that H{t/r) is infinitely differentiable with 
respect to t and that the instantaneous ground state re- 
mains gapped for all time. Following Ref. 1631 we per- 
form successive unitary transformations £A n ) = U^ n \t), 
each attempting to approach with increasing accuracy 
the exact time evolution operator U = U(t); U{— 00) = I. 
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Each transformation is determined from quasistatic con- 
siderations. To define the zeroth-order transformation, 
we first require that it diagonalizes the Hamiltonian at 
each instant of time, i.e., U^HU^ = E (0 K This de- 
termines up to time-dependent phases that we fix 
as follows. Note that propagates the instantaneous 
(zeroth-order adiabatic) eigenstates of H, i.e., 

\i/>P(t)) = U {0) (t)\^°\-oo)). Hence, the additional 
requirement lm.^ip^\d t '4'^) = for all i, which corre- 
sponds to parallel transport,^ determines uniquely. 
Substituting \ip) = \ip (1) ) in idt\ip) = H\tp), we ob- 
tain 



i0 t |^ (1) > = # (1) |V (1) >; 



(10) 



where s = t/r. The term AW = -iU^dgU^, which is 
called the nonadiabatic coupling, is responsible for tran- 
sitions between adiabatic eigenstates. In the basis of 
zeroth-order adiabatic states, it has the elements^ 



A; 



-i(ipi\d s ipj) 

.{i>i\d s H\il}j) 
Ei — Ea 



(j 



(11) 



where the superscripts have been omitted. Due to the 
parallel transport condition, A is purely off-diagonal. We 
now consider the adiabatic energies offfW. Expanding 
the lowest energy eigenvalue of H^ 1 ' with respect to t _1 , 
we find 



E, 



E. 



(0) 



- 2 E 



A {1) A {1) 



p(o) 



E 



(0) 



0{t 



(12) 



The second term is an 0(t~ 2 ) shift of the instanta- 
neous ground-state energy induced by the nonadiabatic 
coupling. Apart from an overall sign, it is identical 
to the induced inertia term that enters in the effective 
Schrodinger equation for the nuclear wave function in 
the Born-Oppenheimer approximation.^ The first-order 
adiabatic eigenstate corresponding to E^ is 



(0) 



A 



(i) 

i0 



^ p(0) 



E 



(0) 



(0) 



0(t~ 



(13) 



the asymptotic series for the wave function in powers of 
t -1 gives a way to derive systematic corrections to the 
ground-state reconstruction. 

IteratiorP^ of the above diagonalization procedure to- 
gether with the definitions — {j( n ) and 

gives the nth-order Hamiltonian 

fj(n) = ft(n-t) _ iT -xfj{n-i)\ ds fj{n-i) _ xhe se quence 

of unitary transformations can be understood as an at- 
tempt to transform the Schrodinger equation to a basis in 
which the nonadiabatic coupling is as small as possible. 
The nth-order adiabatic state | ipo^ ) is an approximation 
for the wave function (up to an overall phase) with er- 
ror of order r~( n+1 ) for all time. However, as described 

in Ref. I63L the sequence of approximants | ipQ 1 " 1 ) cannot 
converge uniformly to the exact solution (it ultimately di- 
verges as nl), for if it did there could be no nonadiabatic 
(Landau-Zener-type) transitions since every member of 
the sequence is asymptotic to the same zeroth-order adi- 
abatic eigenstate at t — oo as at t = — ooP^The approxi- 
mant | ipQ 1 ' ) for which the error is minimum may provide 
an accurate approximation for the wave function, but 
it does not describe nonadiabatic transitions. Assuming 
that the time dependence of H is sufficiently smooth, 
nonadiabatic transitions are nonperturbative, i.e., be- 
yond any power of t _1 . 

Let us consider the IONR approximation from the per- 
spective of the above asymptotic analysis. The lowest 
eigenvalues of the sequence H^ n > provide a sequence of 



adiabatic "ground-state" energies E { n) = E£-'(t). The 
zeroth-order energy Eq ^ — 25q (t) is simply the instan- 
taneous minimum of E v [^] for v — v(t). In a forthcom- 
ing article,^ we shall show that under certain conditions 
there is a sequence of energy functional £$^[7] corre- 
sponding to the Eq 1 \ In contrast to the HK-like energy 
functional E v [7] , for which we can subtract away the con- 
tribution of the external potential energy ^[7], leaving a 
universal functional F[j] — E v [y] — V\y] that does not 
depend on v, the higher-order energy functionals have 
a nonlinear dependence on v and its derivatives with re- 
spect to time, so that the the contribution of the external 
potential cannot be separated. For large enough r, each 



Ey 1 ^ [7] satisfies a local minimum principle yielding Eq 
for 7 = 7''"), where 7'-™^ is the one-matrix correspond- 
ing to IV'o™^)- The 7*-™) are approximations with error of 
order r~( n+1 ) to the exact time-dependent one-matrix. 

The^ n) [ 7 ] can be thought of as defining adiabatic en- 
ergy surfaces. In principle, the constrained minimiza- 
tion of the IONR approximation could be performed on 
one of these higher-order energy surfaces. Since nonadi- 
abatic transitions are not captured in the above asymp- 
totic analysis, we might expect them to be missing in 
the IONR approximation as well. However, the effects 
of nonadiabatic transitions are partially accounted for 
through the effective Schrodinger equation for the nat- 
ural orbitals, and the results that we shall present in 



(n) 



In the absence of magnetic fields, H is real and all of the 
zeroth-order eigenstates i/j^ can be chosen to be real. 
Therefore, ip^ h as the same density, to order r _1 , as 
the state ip^ because the nonadiabatic coupling A 1 - 1 is 
purely imaginary. Moreover, since tp^ is assumed to 
be nondegenerate for all time, it must have everyw here 
vanishing current. The terms containing A 1 -* in Eq. (13) 

generate a current of order t _1 . If tp^ were used to 
construct the two-matrix, Eq. ^ would generate time- 
dependent occupation numbers. In principle, developing Sec. Ill C suggest that nonadiabatic transitions are, in 
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fact, fairly well represented in the IONR approximation. Hermitian matrix Acji = (</>j|A</>i), Eq. (16 1 becomes 



E. Error estimate for the IONR occupation 
numbers 



We can use the hrst-order energy surface to estimate 
the error in the IONR occupation numbers. Let the one- 
matrix obtained in the IONR approximation be denoted 
7. To estimate the error in the occupation numbers we 
take 7^ as a reference and calculate the linear deviations 
with respect to it. Since the error in 7W is 0(t~ 2 ), the 

absolute error of the fi is \fi - f^ I + 0(t- 2 ). We define 
the functional [7] = £^[7] - J dx[y(x, x) - N], 
where fjP'' is a Lagrange multiplier that maintains the 
total particle number N. We assume g4 X) [7], at its 



minimum, satisfies the stationary condition SGv 



(i) 







for an arbitrary variation #7, apart from variations in 
occupation numbers that are exactly or 1. Natu- 
ral orbitals with occupation numbers exactly or 1 are 
called pinned states, and the energy need not be station- 
ary with respect to the variations of pinned occupation 
numbers ! 61 ^ 65 ! 66 ! Under the above assumption, we have 
the stationary condition (for all t) 



dG y v 



(i) 



df k 



(14) 



7=7 (D(t) 



for all unpinned occupation numbers. Similarly, in the 
IONR approximation, we have the stationary condition 



dG v 



dfi 



0, 



(15) 



7=7(i) 



where G v = E v — /1 J dx[ r y(x, x) — N]. Shifting the eval- 
uation point in Eq. ( Jl4| from 7 < - 1 '(t) to j(t), we obtain 
to lowest order 



= T 



dfk 



E 



d 2 E v 1) 
dfkdfi 



(*)-/<(*) 



•£/<** 

i 



(i) 



d SE y v 
df k 5(j>i{x) 



(!) r 



d 6_K 



4>f\x,t)-4>*(^t) > (16) 



where f^ 2 - 1 = t 2 (G v — G v ) and all derivatives are eval- 
uated for 7 = 7. Furthermore, in all second derivative 
terms we have been able to replace Gv by because 
the shift A7 = 7 — 7W is number conserving. Defining 



A/i = fi—f^, A(j>i 



and introducing the anti- 



dfk 



7=7(4) 



^Xkkji fi 



(17) 



We have used 
S 2 E^ 



5j{x' 1 x 1 )Sj(x 2 x 2 ) 



5 2 E V 



d'y(x' 1 xi)Sj(x' 2 x 2 ) 
-X~ 1 (xix / 1 x 2 x' 2 ), 



(18) 

where x _1 is the inverse static response function of the 
instantaneous ground state. Here, Ev has been re- 
placed by E v and the evaluation point has been changed 
from 7 to 7^ - 1 , which both introduce higher-order errors. 
In the basis of natural orbitals, the response function 
x(xix'i, x 2 x' 2 ) = 8^{x\x' 1 ) / 8v{x2X 2 ) is expressed as 



Xijki = J dx l dx l 1 dx 2 dx' 2 cf)*(x 1 )^ j (x[) 

X x{xix' x ,X 2 x' 2 ) <pk{x' 2 )(j)*{X2). 



(19) 



Equation (17) is a set of linear equations relating the 
linear deviations A/j to Acy. It is coupled to another 
set of linear equations obtained from an analysis of the 
linear deviations in the effective Schrodinger equation. 
Together these constitute a linear system of equations 
that can be solved for the A/j and A0j. 

An important simplification can be obtained by con- 
sidering the case that the electron-electron interaction 
is weak. Suppose a coupling constant U is introduced 
into the Coulomb interaction. Consider the following 
four subblocks of X' 1) subblock AA, Xiijji 2) subblock 
AB, xukl (k ^ I), 3) subblock BA and Xklii (k 7^ 0> 
and 3) subblock BB, Xijki {i 7^ 3, k ^ I). If the sys- 
tem remains gapped in the limit U — > 0, subblock AA is 
C(C/ 2 ), subblocks AB and BA are ©(JJ 1 ), and subblock 
BB is 0(1). Therefore, Xnjj = C(C/- 2 ), while Xull and 
Xkm are ©([/- 1 ) and X Zm is O(l). When only the for- 
mer are retained, Eq. (17) gives the following estimate 



for the error of the IONR occupation numbers: 

5(^ 1} - E v ) 



k 

0(U 2 /t 2 ) 



df k 



(20) 

where X = (x -1 ) -1 and (x'^ijkl = SijhlXijlv The 
functional E v — E v is a "warping" of the zeroth-order 
energy surface, and X can be interpreted as an effec- 
tive response function for the occupation number degrees 
of freedom. Since the error of 7 « is also 0{U 2 /t 2 ), 
the absolute error in the IONR occupation numbers is 
0(U 2 /t 2 ). For the model system that we consider in 
Sec. |III| we have verified numerically that the error is 
indeed 0{U 2 /t 2 ). 
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III. 



APPLICATION TO A MODEL SYSTEM 



where 



In this section, we apply the IONR approximation 
and the adiabatic extension approximation in TDDFT 
(ADFT) to a simple model system. By varying the Hub- 
bard parameter U and the characteristic time scale r of 
the external potential, we can study the interplay of elec- 
tron interactions and nonadiabatic dynamics. 

We consider the two-site Hubbard model with two 
electrons.^ The Hamiltonian is 



H 



where cJ CT and Ci a are the creation and annihilation oper- 
ators of an electron with spin a in site i, U — U(ni^hi^ + 
n2-ffi2i), and we have used V±/2 instead of the usual no- 
tation —t for the hopping parameter. In this model, the 
analog of the local external potential v(r,t) is a time- 
dependent bias V3 = V 3 {t/r). It will be convenient to 
write the Hamiltonian as 



4 



1 



C\ 



V 3 (nx 



n 2 



U 



(21) 



(22) 



P 2 = 



-9Q = U 2 + W 2 

-MR = U \2U 2 + 9{V 2 - 3V3 2 ) 



and 



R 



2p 3 



(24) 



(25) 



The definitions of the variables Q and R are conventional 
for cubic equations. The instantaneous eigenstates ji/'i) 
can be expressed as 

- [\W\ 2 +E l {U ~V 3 -Ei)] e" 
\i> i )=N i \ \W\(U-V 3 -Ei) 

\W\ 2 e tV0 

" (26) 

where W = ^(Vi - iV 2 ) = \W\ e~ lVo and N t is the nor- 
malization factor. For completeness, we have recorded 
here the eigenenergies and eigenstates for general V. The 
instantaneous eigenenergies for U = 1, V\ = —2 and 
V3 = t are shown in Fig. [3j 



The operator <x is not re- 



where V = (Vi, 0, V3) and a = (01 , &2, (J3) is the second- 
quantization representation of the Pauli matrices, i.e., 

ft = Eo-(4t> cL) 0- * ( f 2 ° a 

lated to physical spin; below, it will be identified with a 
Bloch pseudospin. In this paper, we consider only con- 
stant Vi and V2 = 0. Generalizing Eq. (22 1 by letting V\ 
and V2 be time-dependent functions is roughly analogous 
to introducing a time-dependent vector potential in the 
case of continuous variables. 

We assume that the initial state is a spin-singlet 
state with S z = 0. As the external potential is 
spin-independent and the singlet and triplet sectors 
are decoupled, this spin configuration will be preserved 
for all times. In the basis {cj^ej^ |0) , 



4-1-4.1.) ^ ' 4f4+ 1°) }> the Hamilt oman is 



V2 [ n c 2; 



U + V 3 ^(Vi-iV 2 ) 

H = I 72^ + lV ^ M Vl ~ %V2) 

jt(Vi + iV 2 ) u-v 3 

The instantaneous eigenenergies, which are also referred 
to as adiabatic energy levels, are the following roots of 
the cubic secular equation: 



E x = 
E 2 = 
E 3 = 



U-pcosQ 
U — p cos — 



U — p cos 





2tt 


3 




u> 


47T 


3 





(23) 



A. Instantaneous occupation number relaxation 
approximation 

The wave function of any two-electron spin-singlet 
state can be factored into a symmetric spatial function 
and a singlet spin function. Thus, it is sufficient to con- 
sider the spatial one-matrix (hereafter, just one-matrix) 
defined as 



(27) 



In the present model, the one-matrix is a Hermitian 2x2 
matrix that can be represented by a Bloch pseudospin 
vector 7 = (71,72,73) according to 7 = I + 7 • a. The 
pseudospin vector 7 should not be confused with the 
pseudospin vector of a two- level system, because I7I is 
time dependent and generally different than 1, while the 
modulus of the latter is always equal to 1, i.e., it remains 
on the Bloch sphere. The natural orbitals in the site basis 
are expressed as 



cos(6»/2)e-^/ 2 
sin(6»/2)e^/ 2 

~sm(e/2)e- ilfi / 2 
cos(6»/2)e^/ 2 



(28) 



and the occupation numbers are f a = 1 + A and /(, = 
1 — A with A = |t*|. In spherical coordinates, 7 = 
A(sin 8 cos ip , sin^siiK/?, cos 9). The 71 component is 
proportional to the kinetic energy, while 72 can be inter- 
preted roughly as an analog of current. Since 73 is local 
in the site basis, it represents the density variable. 



1. Equation of motion 
For the present model, the one-matrix EOM becomes 

id tl = [h n } + u, (29) 

where 7, h and u are 2x2 matrices and h = V ■ <r/2. 
The inhomogeneous term u, which depends on the two- 
matrix (see Eq.[9]), embodies the contribution of electron- 
electron interactions. In the IONR approximation, the 
exact two-matrix functional is approximated by the adi- 
abatic extension of the ground-state functional r[7]. For 
two-electron systems in spin-singlet states, an exact ex- 
pression for the ground-state wave function in terms of 
natural orbitals and occupation numbers is knowrPsl U p 
to sign factors that should be chosen to give the absolute 
minimum of the energyP^ In the present model, |\&) = 



|$ aa ) +v^/h/2 |*w,), where |$«) 



|0) and 



77 is a sign factor that is equal to —1. Here, a\ a and a ia 
are the creation and annihilation operators for natural 
spin-orbital 4>i a (i = a.b). Therefore, the exact ground- 
state reconstruction r[7] = is known and can be 
used to approximate u. In the natural orbital basis, the 
off-diagonal elements of u are found to be 



Uab 



£7(1 + cos/3)sin( 



cost 



(30) 



where f3 — sin A. To calculation u a 6, it is conve- 
nient to use the following expression, which follows from 
Eq. (9): Uij = 2^ a m [a^h^U] |*)- The factor 



cos P in Eq. ( 30 ) represents the occupation number de- 
pendence, while 9 and (p represent the dependence on 
the natural orbitals. In the adiabatic extension approxi- 
mation, the diagonal elements of u are identically zero 
even though we are using the exact ground-state re- 
construction. As a result, the adiabatic extension ap- 
proximation predicts time-independent occupation num- 
bers (cf. Eq. [8|. This deficiency is corrected below in 
the IONR approximation. As an aside, it is interest- 
ing to examine how time-dependent occupation numbers 
are generated in the exact equation of motion. The ex- 
act time-dependent wave function can be expressed as 
I*) = e^' 2 ^jj2 \<& aa ) - e-^l 2 ^fhf2 |* 66 ), where the 
only difference with respect to the ground-state wave 
function is the relative phase factor between the terms. 
If this expression is used to calculate the exact u, one 
obtains u a & = — u'l a = £7(1 + cos /3) sin cos and 
Uaa = —Ubt, = iU cos (3 sin 2 9 sin^>. The exact u generates 
time-dependent occupation numbers because its diagonal 
components are nonzero when there is a nontrivial rela- 
tive phase between the different configurations (Slater de- 
terminants) that comprise the wave function. However, 
no functional approximations are known for the relative 
phases. 

The following EOM for the pseudospin vector 7, which 
follows from Eq. ( 29 ) , provides a geometric interpretation 




FIG. 1: Bloch sphere and trajectory of the pseudospin 7. 
Linear-time potential V = (— 2,0,3£), U — |, t (E (—00, 4.7). 
The point of intersection shows 7 at the time t w —0.261, 
where I7I is minimum; the inner sphere shows |7|min ~ 0.633. 



adiabatic extension approximation: 
9j7 = V x 7 + £7, 



(31) 



where £7 is defined by u — £7 ■ a. Equation (31) is simi 



lar to the Landau-Lifshitz-Gilbert equation or the Bloch 
equation with dissipation. The inhomogeneous term £7 
is responsible for changing the modulus of 7 (Fig. [I]), 
which corresponds to changing the occupation numbers. 
However, in the adiabatic extension approximation, £7 is 



always perpendicular to 7 so that Eq. (31 1 preserves I7 



2. Effective Schrddinger equation 

In the IONR approximation, A(t) is determined from 
the instantaneous minimization of the ground-state en- 
ergy functional E v [7] subject to the constraint that 9 and 
ip are equal to 9(t) and tp{t), the values obtained from the 
effective Schrodinger equation, Eq. ([7]). For the present 
model, 

E v [7] = V ■ 7 + U sin 2 ^ + £7 cos 2 ^ cos 2 9, (32) 

and the value of (3 = sin -1 A that minimizes the energy 
for given 9 and ip is 



cot 



£7 sin 2 9 
2 \f . a 



(33) 



The effective Schrodinger equation for the orbitals, 
Eq. Q, becomes 



for the conservation of the occupation numbers in the 



id t (f>t = (i + V e s) <j>i, 



(34) 
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where t = V\<j\j2 and v c g = v + v ee with v = V3CT3/2 
and v ee — V ee ■ a/2. In spherical coordinates, V ee has 
components 



V ee ,e — ~U cot — sin 9 cos 9 

= 0. 



(35) 



The radial component of V ee (the component parallel 
to 7), which corresponds to the "diagonal" component 

We 



v e e,bb, is nonunique as discussed in Sec. II B 



remark in passing that v ee , apart from this nonunique 
component, is equal to SW^/Sj, where W[j] is the 
ground-state electron-electron interaction energy func- 
tional. Therefore, in the IONR approximation, v ee is 
the adiabatic extension of the electron-electron part of 
the generalized ground-state KS potential.^ 



3. Linearization 

It is interesting to examine the role of the nonlinear- 
ity of V ee [7] by carrying out a linearization and semilin- 
earization with respect to a reference dynamics or zeroth- 
order dynamics. A suitable reference is the instantaneous 
ground-state one-matrix J^°'(t), from which the exact 
one-matrix does not deviate too greatly in the adiabatic 
regime. Thus, the lo west -order (linear) approximation 
consists in solving Eq. (34 1 with V ee evaluated at 9 = 9^ 



and tp = tp(°\ where 9^ v> and tp 1 -- ) are the angular vari- 
ables obtained from 7^°' (<^ ' = if V2 = 0). In the 
next-lowest order (semilinear) approximation, we solve 
Eq. ( 34 ) self-consistently with the potential 



ysl _ ycffj 7 (0)j + 



dV 



off 



dcos9 



5 (cos 9) 



dV 



cos tp 

where, for example, 

(5(cos 9) = cos 9 — cos 

0(0) 

= (M 2 -M 2 ) 



8 (cos tp) 



(|<4 0) | 2 -|4T), (37) 



(°)|2^ 



where are the elements of <j) a . We have linearized 
in the variables cos 9 and cos tp rather than 9 and tp 
because cos 9 is closely related to the density variable 
73 = A cos 9, which, in the following section, will be used 
in the linearization of the TD KS equations. The linear 
and semilinear approximations are compared with the 
full IONR approximation in Sec. |III C 5| 



B. Adiabatic extension in time-dependent density 
functional theory 

To better understand the nature of the IONR approx- 
imation, it is helpful to compare it with ADFT. The 



methods are similar because in both cases the orbitals 
are determined by an effective Schrodinger equation, in 
which the effective potential is the adiabatic extension of 
a ground-state effective potential. 



1. Time dependent Kohn-Sham equations 

For the two-site Hubbard model, the TD KS equations 
in the adiabatic extension approximation are 



id t <t> a = (i + v ae )^ a , 
73 = 2 I -<7 3 |<?!>a>, 



(38) 



where the factor of 2 comes from the fact that the orbital 
(p a is doubly-occupied and v ae = (V ae /2)a 3 with 



V ae (t) = V(t) + V H (t) 

= _ V1 •»(«) 



dE xc 
973 



73=73(4) 



Fx cot [6 s (t)}, 



(39) 



where Vh (t) is the time-dependent Hartree potential and 
E xc is the ground-state xc energy functional. The adia- 
batic KS energies are e a ,b = z f(Vi/2) csc9 s . The density 
variable in this model is 73 = cos 9 S , where 9 S should not 
be confused with 9 defined in Sec. IIII Al 



2. Linearization 



In analogy with Sec. Ill A 3 we carry out a lineariza- 
tion and semilincarization with respect to the instan- 
taneous ground-state density 73 (*)■ The lowest-order 



(linear) approximation consists in evolving Eq. ( 38 ) with 



(36) V£ 



(0) 



In the next-lowest order (semilinear) 



approximation, we evolve Eq. ( 38 ) with 



K S i(73)=Ke(7f ) + 



^73 



(0) 
7 3 



<$73 



= K e (7 3 U, ) + (x7 1 -X- 1 )^73 (40) 

where Xs and x are the instantaneous static ground-state 
response functions. Thus, the better Xs approximates Xi 
the less significant is the nonlinearity. The linear devi- 



ation of the density is 673 = | (ai | 2 — \a 2 \ 



,(°)|2 



(K 



(0)|2 



| 2 ), where Oj are the elements of the KS orbital 



The TD KS equations become 



id. 



Oi 
(to 



C^73)o- 3 ] 



a 1 

«2 



(41) 



where C = xj 1 — X • Although the potential has been 
linearized with respect to ^73, the equations remain non- 
linear in the amplitudes a.j. This nonlinearity is essen- 
tial for the accurate description of nonadiabatic effects. 
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FIG. 2: (Color online) Time-dependence of the one-matrix. Linear potential, r = 1, U = 1. Top left, top right, and bottom 
left: 71, 72, and 73. Bottom right: adiabatic many-body energies [solid lines] and adiabatic KS energies [dotted lines]. 



Equation (41 1 has the form of a nonlinear LZ problem. 



In fact, similar equations have been studied as simple 
models of the Gross-Pitaevskii equation in the context of 
trapped Bose-Einstein condensatesP^HH] The linear and 
semilinear approximations are compared with the full 
ADFT in Sec.lIIIC5J 



C. Numerical simulations 

We have performed simulations comparing the IONR 
approximation, TDHF, ADFT and the numerically ex- 
act solution. The following two types of time-dependent 
external potentials were considered: a linear potential 
V3 = t/r and a pulse potential V3 = 5/ cosh(£/r). The 
case of the linear potential can be interpreted as the scat- 
tering of two species — one binding two electrons and one 
that can accept electrons. We are interested in the dy- 
namics of systems that start in the ground state, and in 
principle we can take any time as the initial time to . Since 
for both potentials t — ~ 00 is a point where the nonadia- 
batic coupling vanishes and the instantaneous eigenener- 
gies are nondegenerate, it is intuitively clear that ~/(t) will 
approach a unique time-dependent function as we let to 
tend to —00, always choosing the ground state for the ini- 
tial state. In our simulations, we set to = —p and choose 
larger and larger values of p until "f(t) is sufficiently con- 
verged for all t > to. We discuss results obtained for the 
linear potential in Sees. |III C lpll C 5| and for the pulse 
potential in Sec. Ill C 6 



1. Nonadiabatic effects 

The time-dependence of the one-matrix is plotted in 
Fig. [2] for U = 1, Vi/2 = -1 and V 3 = t (r = 1). For 
these parameters, the system is in an adiabatic regime 
and, as seen in the plot of 73, both electrons are trans- 
ferred smoothly from site 1 to site 2 (73 = 1 — > —1). The 
model has two dimensionless parameters: U /V\ and V\T. 
In all of the following, we take the hopping parameter 
Vi/2 to be equal to —1, and refer to U and r as dimen- 
sionless parameters. Figure [3] shows the deviation of the 
one-matrix in Fig. [2] from the instantaneous ground-state 
one-matrix 7^(t). We can identify two principal nona- 
diabatic effects: 1) coherent oscillations in j(t) as t — > 00 
and 2) mixing with excited states near avoided crossings 
of the adiabatic energy levels. The asymptotic oscilla- 
tions originate from nonadiabatic (LZ) transitions. These 
are transitions between the initial and final (t — > ±00) 
adiabatic eigenstates and generally leave the system in 
a nonstationary state. The admixture of excited states 
near avoided crossings is apparent in the deviation of 72 
from 7| 0) ( 7 f ) = for all time) and A from . The 
projection of the wave function onto an excited state can 
be quite large near an avoided crossing, even if ultimately, 
for large times, the amplitude becomes very small. The 
right panels in Fig. [3] show that the instantaneous corre- 
lation energy, E c = E — E^f 1 is strongly correlated with 
the deviation of the occupation numbers from or 1. 

The interplay between electron-electron interactions 
and nonadiabatic effects is, in leading order, mediated 
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FIG. 3: (Color online) Time dependence of the one-matrix relative to the instantaneous ground-state one-matrix. Linear 
potential, r = 1, U = 1. Top left and bottom left: 71 — 7^°' and 73 — Top right: A. Bottom right: correlation energy E c . 



by the adiabatic energy levels, which will be discussed in 
detail below. Figure [4] shows the one-matrix in a strongly 
interacting case; r = 2 and U = 4. Stronger interactions 
lead to more pronounced curvature in the energy level 
profiles (compare Figs. [2] and [4|. Interestingly, the adia- 
batic KS and HF energies do not display such curvature. 

In Fig. [l] we have plotted the trajectory of the pseu- 
dospin vector 7 for U = 7/2 and r = 1/3. The vector 
starts at the north pole at t = —00 and, in following the 
driving vector —V, it rotates toward the south pole. Due 
to the finite probability of nonadiabatic transition, it does 
not reach the south pole and instead spirals around it per- 
petually. Near the equatorial plane of the sphere (t ~ 0), 
the occupation numbers deviate significantly from their 
original values [f a (— 00) = 2 and /&(— 00) = 0] so that 
7 leaves the surface of the Bloch sphere (|7|=1). Re- 
markably, the trajectory exhibits a fairly regular spiraling 
pattern even at such times. In contrast to the persistent 
spiraling as t — > 00, for times near the avoided crossings 
the spiraling is not tangential to the surface of the Bloch 
sphere. 



2. Asymptotic oscillations 

By studying the amplitude and phase of the asymp- 
totic oscillations, we can extract information about the 
nonadiabatic transitions. For the range of parameters 
investigated here, the asymptotic oscillations are dom- 
inated by a single time-dependent frequency f2(t) = 



E 2 {t) — E\{t). Thus, the one-matrix is well described 
by the expression 



7i = 7,; + A, cos (j dt'VL{t') - 9, 



(42) 



The quantities 7 y i , Aj and Qi are related to the final am- 
plitudes of the adiabatic states. Let us write the many- 
body wave function as 

3 

l*(*)) = Y,^ t ) e ~ iS ° dt ' Eh{t ' ) \Mt)) , (43) 
k=l 

where \ipk) and Ek are the instantaneous eigenstates and 
eigenenergies given in Eqs. (23 1 and (26). Except for low 



values of r (r < 1/4) , the final amplitude of the highest 
energy adiabatic state ^3) is much less than those of the 
lowest two states because it is separated from the ground 
adiabatic state by a larger energy gap. Therefore, as a 
first approximation in the adiabatic regime, let us trun- 
cate the expansion in Eq. (43) to the lowest two levels. 
Then, we find 

7i = N a ^i|^i|^> + |ca| 2 <^|^i|^) 
1 



Aj = 2|ci||c 2 | <Vx j ifi 1^2) 
&i = Arg(c 2 /ci). 



(44) 



DC. 



where here c%, C2, \ipi) and \ip 2 ) ar e evaluated at t 
In Fig. [5], the amplitude |Ai| and phase 81 obtained from 
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FIG. 4: (Color online) Linear potential, r = 2, U — A. Top left and top right: 71 and 72. Bottom left: A — A^K Bottom 
right: adiabatic energy levels (colors same as Fig. [5] with the addition of HF eigenvalues [dashed lines]). 



fitting the simulations to Eq. ( 42 ) are shown as functions 
which can be interpreted as a scattering "veloc- 



of 



ity." The amplitude of the oscillations is found to be 
quite sensitive to r _1 , even "collapsing" for particular 
values. Similarly, if |Ai| is plotted with respect to U for 
fixed t -1 , it is seen that there is a series of values of U for 



which the oscillations collapse. In Sec. Ill C 3 it will be 



shown that both of these instances of collapse are mani- 
festations of the same phenomenon, namely an interfer- 
ence effect often referred to as Stueckelberg oscillation^ 
or Landau-Zener interferometry. Here, oscillations refers 
to the fact that the effect is quasiperiodic and mediated 
by a sine factor. Stueckelberg oscillations are observable 
in a variety of settings (see for example Refs. 1701 and I7T|) . 
As far as we can ascertain within the numerical preci- 
sion of the simulations, the amplitude of the oscillations 
collapses all the way to zero in all of the approximations 
investigated — TDHF, ADFT and IONR. In contrast, 
in the numerically exact solution the amplitude of the 
oscillations reaches instead a finite minimum value. 

The results in Fig. |5j taken as a whole, indicate that 
the IONR approximation and ADFT perform compara- 
bly. The IONR approximation gives more accurately the 
critical values of r^ 1 for which the oscillations collapse. 
Although in principle TDDFT need only give exactly the 
diagonal element 73 — the density variable in the present 
model — and not the off-diagonal elements 71 and 72 , it 
is nevertheless appropriate to compare the methods on 
the basis of the asymptotic oscillations in 71 because 
these oscillations have the same origin as correspond- 
ing oscillations in 73, which, however, have an ampli- 



tude that decays to zero as t — > 00 due to the fact that 
lim t _ i . 00 (V>i|<73|V'2) — 0- The phase Oi exhibits a strik- 
ing resonance behavior with respect to t _1 , jumping by 
7r as it passes through ±ir/2. The resonances coincide 
with the minima of |Ai|. In the elucidation of this reso- 
nance phenomenon, it will be helpful to have a method 
for estimating the asymptotic final values of c\ and C2 
that appear in Eq. (44). We shall now describe such a 
method. 



3. Independent crossing approximation 

The final amplitudes Cj(oo) of the adiabatic states are 
related to the initial amplitudes c.;(— 00) by a unitary 
scattering matrix. In the adiabatic regime, the scattering 
matrix can often be calculated with good accuracy in the 
so-called independent crossing approximation (ICA).^l 
Consider a multilevel system in which the adiabatic en- 
ergy levels undergo pairwise avoided crossings. In the 
adiabatic limit, nonadiabatic transitions are typically lo- 
calized near the avoided crossings and each avoided cross- 
ing can be described by a scattering matrix that connects 
the amplitudes of its incoming and outgoing adiabatic 
states. Between adjacent avoided crossings, the evolu- 
tion is nearly adiabatic and the components of the wave 
function simply acquire dynamical and geometric phases 
(there is no geometric phase in our model because V is in 
the xz-plane and does not encircle the origin; however, 
it will appear when V in fully 3-dimensional and sub- 
tends a nonzero solid angle). An approximation to the 
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scattering matrix for the full time interval can be con- 
structed by taking the time-ordered product of the indi- 
vidual scattering matrices, interposed by diagonal matri- 
ces describing the dynamical and geometric phases, pro- 
vided such a time-ordering is possible. In the case of 
the linear-time external potential, the ground adiabatic 
state undergoes two avoided crossings with the first ex- 
cited adiabatic state (see Fig. [4j. These crossings are 
related by reflection about t — 0. Although the adia- 
batic energy levels do not intersect for real time, there 
are generally crossing points in the complex time plane. 
For example, for U = 1 and r = 1, levels 1 and 2 cross 
at t rj ±1.2737 ± z2.4584. Let t a be the crossing point 
with Rc t < and Im t > 0, and let tb be the crossing 
point with Re t > and Im t > 0. Between the levels 2 
and 3, there are two complex conjugate crossing points; 
for U = 1 and r = 1, they are t = ±il.0757. Let t c be 
the crossing point in the upper half plane. Levels 1 and 
3 do not cross for any complex time, so direct transitions 
between these levels are expected to be very weak. 

The scattering matrix for the avoided crossing near t a 
can be approximated by 



Ra ~T* a 

S a ~ \ T a R* a 
1 



(45) 



Due to the time reversal symmetry of the adiabatic 
energy levels and the nonadiabatic coupling A 2 \ — 
—i {ip2\dtipi) between levels 1 and 2, the corresponding 
matrix near tb is S b = (S a y . The scattering matrix near 
t c is approximated by 



1 
,S" = | R c —T* 

T c R* 



(46) 



In the ICA, the full scattering matrix connecting the adi- 
abatic states at t = — oo and t — oo is 



S = D ob S b D bc S c D ca S a D 



ca qcl T~\aO 



(47) 



where D a P are diagonal matrices containing dynamical 
and geometric phases, for example, D a0 has the elements 
Df k = e - m ti a dt(B h -i(^ k \d t ^ k )) _ The fml scattering ma _ 

trix gives an estimate for the asymptotic final amplitudes 
4 CA = S n : 

c{ CA (oo) = \R a \ 2 +R c \T a \ 2 e-^ 

c I 2 CA (oo) = -R a T a (e^ 2 ~R c e-^ 2 ) 

4 CA (oo) = T a T c e~ m fo a dt'E le -ifHf: : dt'E 2e -ifHfl dt'E 3 

(48) 

where $ = Re f* b dt Q(t). In the limit t" 1 ->• 0, T c ->• 
and R c -> 1. Thus, for low r _1 , c| c ' A (oo) 0, and 
c[ CA (oo) and c^ CA (oo) can be approximated by a two- 
level (TL) approximation 

cT L (oo) = \R a \ 2 + \T a \ 2 e- 1 * 



TL 



(oo) 



-i2R a T a sin 



4> 




FIG. 5: (Color online) U = 1. Amplitude |Ai| and phase 
Bi of the asymptotic oscillations of 71. Arrows and vertical 
lines represent the resonances. Dashed vertical lines repre- 
sents the reso nances predicted by the two-level approxima- 
tion, Eq. (l49l). 



(49) 



Stueckelberg oscillations are encoded in the sine factor, 
which describes the collapse of the oscillations and the 
accompanying phase jumps. The amplitude of the os- 
cillations I Ai| is proportional to |c2(oo)| and, in the 
two-level ICA, 02(00) vanishes for Re J t b dtCl(t) — 27m; 
n = 1,2,3,... In this Bohr-Sommerfeld-like condition, 
fl(t) = E2(t) — E\(t) depends implicitly on U and r _1 
(and, generally, on any other parameters of the system 
that affect the adiabatic energy levels). Therefore, the 
Stueckelberg oscillations with respect to U have the same 
origin as those with respect to r _1 . The vanishing of 
c^ L (oo) can be interpreted as the destructive interference 
between the two different "pathways" through the energy 
level diagram that start in level 1 at t = — 00 and end in 
level 2 at t — 00. In the first, a nonadiabatic transition 
occurs near t a , in the second, near t b . The phase asso- 
ciated with each pathway is the product of three types 
of phases — dynamical, geometric and scattering. In 
the present model, as already mentioned, the geometric 
phase is zero. And, due to the condition S b — (S a Y , 
the contributions of the scattering phases contained in 
S a and S b cancel out. 

The phase of the oscillations Oi depends on the relative 
phase of ci and c 2 (here and in the following we omit 
the argument t = 00). To the extent that c 2 can be 
approximated by c 2 L , we expect to see perfect step-like 
jumps of 7T in 0i, coinciding with the sign changes of 
sin(<i>/2) as it passes through zero. In Fig. [5] the jumps 
predicted in the two-level approximation are indicated by 
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dashed vertical lines. In contrast to the numerically exact 
result [solid black line] , all of the approximations display 
discontinuous phase jumps. Furthermore, the phase in 
all of the approximations is advanced with respect to the 
exact result. 



4- Effect of the third adiabatic state 

The numerically exact solution in Fig. [5] displays a 
resonance behavior in which the phase increases rapidly 
but smoothly by it at particular values of t -1 . To un- 
derstand the broadening of the resonances (note that it 
is not the nominal amplitude |Ai| but rather the quan- 
tity |Ai| _1 that becomes very large on resonance), it is 
necessary to take into account the third adiabatic state 
l^). The broadening of the resonances is primarily due 
to transitions to the third state, which destroy the per- 
fect destructive interference that occurs in the two-level 
approximation. 

We shall now calculate the width of an arbitrary reso- 
nance and show that it is proportional to the minimum 
value of |Ai|. We begin by expressing c[ CA and c I 2 CA as 
follows: 

c{ CA = l-\T a \ 2 e- i (*-P°y 2 <; 

ci CA = -R a T a e^l\, (50) 

where p c — Arg(i? c ) and 

C = e i l*- p »V 2 - |i? c |e- l( *-^ )/2 (51) 

is a key quantity for describing the resonances. The ar- 
gument of C> 



Arg(C) = tan 1 



1- R c 



(52) 



rises smoothly from — ir/2 to ir/2 near $ — p c = 2wn. 
The width (in terms of r _1 ) of a resonance at r = t' is 
defined to be |d0i/d(T _1 )|~i; T ,. In the following, we shall 
assume that \T a \ <C 1, which is the case for sufficiently 
low t^ 1 . Hence, the argument of c I 1 CA is approximately 
constant, and the r-dependence of Oi comes primarily 
from the argument of c I 2 CA . We find 



d0i dArg(C) , dAvg(R a T a ) 1 dp c 



d(r-!) d(r-!) 



d{T- x ) 2d(r- 1 )' 



(53) 



The first term is 0(t 2 ) because $ ~ r as r _1 — > 0, so 
that the second and third terms, which are 0(1), are 
negligible. 

To express d0i/d(r _1 ) in terms of |Ai|, we first con- 
sider the condition 



= 



d(r-i) ~ (l-\R c \) 2 



Ai\\Rc\ . ,^ ,d(<$>~p c ) 

sm($-p e ) \. J. CJ , (54) 



d(r-i) 



neglecting terms of order d\cz\ 2 /d{T 1 ). Equation (54) 
gives the resonance condition $ — p c = 2im. The phase 
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FIG. 6: (Color online) Linearization and semilinearization of 
the IONR approximation. Linear potential, r = 1, U = 1.25. 



Pc due to scattering with the third level, shifts the res- 
onances with respect to the positions predicted in the 
two-level approximation. It becomes signifi cant for mod- 
erately large r _1 (see Fig. [5]). Using Eqs. (48) and (54) 
in Eq. (53), we obtain 



d0i 



d(r-i) 



1 1 

2\k;\ 



\T C \ 2 ^ dIAilMr- 1 ) 



\Rr. 



sin($ - p c 



(55) 



This establishes that the resonance width is proportional 
to |Ai|. Within the ICA, the minimum value of |Ai| 
on resonance is V2\T a R a \(l - \R C \), which vanishes very 
rapidly as r^ 1 — > because \T a \ —> and \R C \ —> 1. 
Thus, there is an infinite sequence of increasingly sharp 
resonances in the limit r^ 1 0. 

The asymptotic oscillations can be described exactly 
by including the effect of the third adiabatic state that 
was neglected in the two-level approximation. Using 
limt^oo (tpi \&i\ipi) = and lim^oo (tpi \ai\1p3) = 0, we 
find 7 y 1 = and the following asymptotically exact ex- 
pression: 



7i(i) = NM<Vil*ilV>a>coB / dt'ct>(t')-Ut 



7 21 



+ NN H cos ( / dt'(j)(t') + ut-e 32 

(56) 



where Ci and \ipi 
2b 2 r/t and % = 



are evaluated at t — 00, <j) = tjr 



Arg(c i /c i 



The third state induces 



additional oscillations with time-dependent frequency 
E3 — E2 ~ 1 
while if |c3 



0, Eq. (56) reduces to Eq. (42 1, 



-U. If |c 3 

i |ci|, the third state generates beats with 
period 2w/U. The ICA prediction for the final amplitude 

of the third state, c 3 = 531 = e~ 2jrb r , is actually exact. 
This is a special case of a general result for n-level sys- 
tems with linear time dependence, namely, it has been 
proved^ that the ICA is exact for the scattering matrix 
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elements Si„ and S n x, where 1 and n denote the lowest 
and highest energy adiabatic states. 

In the many-body system, the interference of dynam- 
ical phases arises from nonadiabatic transitions at two 
distinct and well separated times (t a and tb). Therefore, 
it is interesting to ask how ADFT and the IONR approx- 
imation, which are devoid of memory dependence, are 
able to describe such an effect. For ADFT the reason is 
clear: the same interference phenomenon operates in the 
KS system. How nonadiabatic transitions and the inter- 
ference phenomenon are captured in the generalized KS 
system represented by Eqs. §f§ and ([8| is an interesting 
question for future study. Both ADFT and the IONR ap- 
proximation greatly underestimate the resonance width; 
in fact, within the numerical accuracy of our simulations, 
we were unable to resolve any of the resonance widths. 
This deficiency of the approximations is due to the ab- 
sence of a spectator KS state (the KS and generalized 
KS systems have only two states), analogous to the third 
state in the many-body system, that can "scatter" with 
the states participating in the interference phenomenon. 

It is difficult to calculate the scattering matrix el- 
ements in the general case. In the adiabatic regime 
the square moduli of the off-diagonal elements, which 
give transition probabilities, can be estimated with the 
Dykhne formulap^^ e.g., for the transition near t a , 



-2Qf* a dt{E 2 -Ex) 



(57) 



The Dykhne formula is universal in the sense that it de- 
pends only on the adiabatic energies and not the non- 
adiabatic coupling. For the "reflection" component of 
S a , we have \R a \ = \Jl — \T a \ 2 . However, having only 
the moduli of T a and R a is not sufficient to fully describe 
the asymptotic oscillations: the phases are also impor- 
tant. In Eqs. (48), we see that the arguments of R a , 
T a and R c directly influence 8i and the location of the 
resonances. The U and r dependence of these arguments 
is manifest in the nonvanishing slopes of the plateaus in 
Fig. [5] 

Progress in the context of time-dependent functional 
approximations might require confronting memory de- 
pendence in an explicit way. From this point of view, it 
is noteworthy that the Dykhne formula, as well as the 
dynamical and geometric phases in D a ^ , contain a type 
of memory dependence due to the time integrals. The 
functional dependence enters through the adiabatic en- 
ergies, which can be taken to be functionals of the instan- 
taneous ground-state one-matrix (ground state density) 
rather than the true time-dependent one-matrix (time- 
dependent density). However, as mentioned above, the 
Dykhne formula gives only the probability of nonadia- 
batic transition at an avoided crossing, and the scattering 
phases are also important. The nonadiabatic coupling is 
expected to be important in the calculation of scatter- 
ing phases. While the calculation of excited state ener- 
gies within linear response theory is fairly well developed, 
much less is known about the nonadiabatic couplings (for 
example, see Refs. [771 and 175)) . 




FIG. 7: (Color online) Linearization and semilinearization of 
ADFT. Linear potential, r = 1, U = 1.25. 



5. Validity of linearization and semilinearization 

The time-dependent Kohn-Sham potential and the ef- 
fective potential in Eq. Q are nonlinear functionals 
of the time-dependent density and time-dependent one- 
matrix, respectively. By considering the linear and semi- 
linear versions of these equations, we can investigate the 
importance of nonlinearity in the description of nonadi- 
abatic effects. 

The linear and semilinear versions of Eq. ([7]) in the 
IONR approximation were introduced in Sec. |III A 3| 
Figure [6] compares the linear and semilinear versions with 
the full (unmodified) equations for representative values 
of U and r. Even the linear version generates asymptotic 
oscillations, although their amplitude, phase and period 
are generally incorrect. The period is corrected in the 
semilinear and full versions. The semilinear version im- 
proves significantly also the amplitude, phase and mean 
value of the oscillations with respect to the linear version. 



In Sec. |IIIB 2[ similar linearizations were performed 
in ADFT. Figure [7] compares the linear and semilinear 
versions with the full equations. As found in the IONR 
approximation, the amplitude and phase of the asymp- 
totic oscillations are quantitatively incorrect in the lin- 
ear version and the nonlinearity of the semilinear version 
brings them into better agreement with the full version. 
Surprisingly, only this lowest-order nonlinearity is suffi- 
cient to bring the semilinear version into nearly perfect 
agreement with the full version, except for large values 
of U and low values of r (U > 8 and r < |). In con- 
trast to the linear version of Eq. Q, the linear version of 
the KS equations achieves the correct period because the 
difference in the adiabatic KS energies tb ~ £a equals the 
difference E% — E\ in the limit t — > oo (see Fig. [4]). How- 
ever, this is a special property of the linear-time potential 
due to its divergence in the limit t — > oo. 
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FIG. 8: (Color online). Time dependence of the one-matrix for a pulse potential V3 = 1/coshi; U = 3. 



6. Pulse potential 



IV. 



CONCLUSIONS 



Up to this point, all of the simulations we have re- 
ported have used an external potential with linear time 
dependence. As such a potential diverges in the limit 
\t\ — > 00, it has the special property that the ground 
state at t = —00, in which both electrons occupy site 1, 
is an uncorrelated single Slater determinant. To examine 
the performance of the approximations in systems with 
correlated initial states, we have considered additional 
time-dependent external potentials. Here we report sim- 
ulations for a pulse-shaped potential, V3 = l/cosh(i/r). 

In Figure |8j the IONR approximation and ADFT are 
compared with the numerically exact solution. The per- 
formance of ADFT has worsened for this external po- 
tential. The most significant shortcoming of ADFT is 
that the period of the asymptotic oscillations is too short. 
This is the opposite of what one would expect on the ba- 
sis of the adiabatic KS energies: £}, — e a is smaller than 
E2 — E\, suggesting that the period of the KS oscilla- 
tions should be too long. This counterintuitive behavior 
is accounted for by the nonlinearity of the KS equations. 
Although ADFT gives a qualitatively incorrect descrip- 
tion of 71 , this should not be viewed as a deficiency of the 
approximation, as TDDFT is only guaranteed to repro- 
duce the density variable 73 and not the full one-matrix. 

The IONR approximation performs better than ADFT 
for this external potential. Its most pronounced defi- 
ciency is that it does not capture the asymptotic oscilla- 
tions in 71 (i) and A(t). 



Describing strongly-driven electron dynamics from 
first principles is a challenging problem. Simulations of 
a multi-configurational wave function with enough terms 
to adequately describe correlation have so far been lim- 
ited to small systems or short simulation times. Time- 
dependent density-functional theory is a less computa- 
tionally demanding approach capable of treating much 
larger systems. However, approximations for the xc po- 
tential, a complicated nonlocal and memory-dependent 
functional of the density, must be introduced. Essentially 
all calculations to date have employed the adiabatic ex- 
tension approximation. Its known deficiencies in linear 
response calculations may have counterparts in real-time 
simulations. Steps toward overcoming the deficiencies of 
the adiabatic approximation in linear response TDDFT 
has been made^^U by using the one-matrix instead of 
the density as basic variable. 

Working with the one-matrix may have advantages for 
real-time dynamics as well. However, applying the adi- 
abatic extension approximation to the one-matrix EOM 
yields time-independent occupation numbers when any of 
the available ground-state functionals are used. We have 
proposed a simple modification of the adiabatic extension 
approximation, here called the IONR approximation, in 
which time-dependent occupation numbers are obtained 
"on-the-fly" from a condition of instantaneous relaxation 
to the minimum of an adiabatic energy surface. The 
IONR approximation has the advantage that it yields 
time-dependent occupation numbers even for the avail- 
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able functional, which have proved successful for many 
ground-state properties. The motivation for the instanta- 
neous relaxation condition can be understood in light of a 
sequence of adiabatic energy surfaces that arises from an 
asymptotic analysis of the many-body Schrodinger equa- 
tion in the limit r — > oo. Each energy surface obeys an 
instantaneous minimum principle giving an approxima- 
tion to the exact one-matrix with error of order r~( n+1 \ 
However, arbitrary accuracy cannot be achieved on the 
basis of this sequence because the existence of nonadia- 
batic transitions causes the sequence to divergence. 

We performed simulations for a model system, which 
demonstrated that the IONR approximation captures 
fairly well nonadiabatic effects, such as LZ-type transi- 
tions, even though it lacks memory dependence. Non- 
adiabatic transitions leave the system in a nonstation- 
ary state, resulting in persistent oscillations in the ob- 
servables. The amplitude and phase of the oscillations 
undergo resonances when either the electron-electron in- 
teraction strength U or the characteristic time scale r 
is swept through critical values. We have found that 
the IONR approximation describes the persistent oscilla- 
tions qualitatively correctly over a wide range of U and 



r. It becomes quantitively correct for low values of U 
and large values of r. This is remarkable because the 
resonance behavior is due to an ultra-nonlocal interfer- 
ence of dynamical phases and scattering phase shifts in 
the interacting system. 

Although the IONR approximation is able to describe 
the lowest-order nonadiabatic effects, it has an important 
limitation: being based solely on the adiabatic extension 
of ground-state functionals, it does not apply in situa- 
tions where the time-dependent wave function deviates 
greatly from the instantaneous ground-state wave func- 
tion. Thus, while it may be relevant for adiabatic quan- 
tum control problems, it is not applicable to the type of 
quantum control problem studied in Ref. 1481 where the 
target state is an excited state. 
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negation of a Greek negation, nonadiabatic is used in quan- 
tum dynamics to describe corrections to purely adiabatic 
dynamics, i.e., corrections to the dynamics obtained by ne- 
glecting transitions between the instantaneous (adiabatic) 
eigenstates. 

When the Hamiltonian is real and its ground state is non- 
degenerate, the exact reconstruction must be real. 
In the Born-Oppenheimer approximation, the adiabatic 
electronic eigenstates depend implicitly on the nu- 
clear trajectories ~R, n {t), so that r _1 Aij contains the con- 
tribution (ipi\V„ ■ Pn\ipj), where V n = dR n /dt and 
P n = — iV n - In the literature, the factor (V'ilVn^'j) i s 
often referred to as the nonadiabatic coupling vector. 
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